## Callin Switzer ## 29 Nov 2016 ## 6 Dec 2016 Update -- re-digitized pollen and anthers very carefully ## may not need any smoothing ## 8 Dec Update: used cross validation for decide smoothing parameter for smoothing spline ## 8 Feb 2017 Update: Cleaned up code for submission as supplemental info ## 22 June 2017 Update: Made example plots for position/speed/accel ## Kalmia pollen and anther kinematics ### 1. Read in digitized files ### 2. Smooth digitized points, using smoothing spline with tuning parameter from cross validation ### 3. Use smoothed points to calculate velocity and acceleration (normal and tangential)
ipak <- function(pkg){
new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)
}
packages <- c("ggplot2", "scales", "multcomp", "plyr", "car", "lme4", "signal")
ipak(packages)
## ggplot2 scales multcomp plyr car lme4 signal
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE
# print session info and time the code was run last
Sys.time()
## [1] "2017-06-22 20:11:20 EDT"
sessionInfo()
## R version 3.3.2 (2016-10-31)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X El Capitan 10.11.6
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] signal_0.7-6 lme4_1.1-12 Matrix_1.2-8 car_2.1-4
## [5] plyr_1.8.4 multcomp_1.4-6 TH.data_1.0-8 MASS_7.3-45
## [9] survival_2.40-1 mvtnorm_1.0-5 scales_0.4.1 ggplot2_2.2.1
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.9 nloptr_1.0.4 tools_3.3.2
## [4] digest_0.6.12 evaluate_0.10 tibble_1.2
## [7] gtable_0.2.0 nlme_3.1-130 lattice_0.20-34
## [10] mgcv_1.8-16 yaml_2.1.14 parallel_3.3.2
## [13] SparseM_1.74 stringr_1.1.0 knitr_1.15.1
## [16] MatrixModels_0.4-1 rprojroot_1.2 grid_3.3.2
## [19] nnet_7.3-12 rmarkdown_1.3 minqa_1.2.4
## [22] magrittr_1.5 backports_1.0.5 codetools_0.2-15
## [25] htmltools_0.3.5 splines_3.3.2 assertthat_0.1
## [28] pbkrtest_0.4-6 colorspace_1.3-2 quantreg_5.29
## [31] sandwich_2.3-4 stringi_1.1.2 lazyeval_0.2.0
## [34] munsell_0.4.3 zoo_1.7-14
# set directory that contains datasets
dataDirectory = '/Users/callinswitzer/Dropbox/ExperSummer2016/Kalmia/Manuscript/DatasetsSupplemental/'
# read in metadata
dfile <- "LaurelsOnly.csv"
metDat <- read.csv(paste0(dataDirectory, dfile))
metDat <- metDat[metDat$redigitizedFile!= "", ]
# set constants:
fps <- 5000 # frames per second
Note: The smoothing parameter was chosen with 10-fold CV
# read in each .csv file for analysis
digdirect <- "/Users/callinswitzer/Dropbox/ExperSummer2016/Kalmia/Manuscript/DatasetsSupplemental/CleanKalmiaDigitized/"
newDF <- data.frame()
for(ii in 1:nrow(metDat)){
ddfile <- paste0(digdirect, metDat$redigitizedFile[ii])
dp <- read.csv(ddfile)
# calibrate locations, based on digitized pin or other object
# calibration points
pin <- data.frame(dp$pt1_cam1_X, dp$pt1_cam1_Y, dp$pt2_cam1_X, dp$pt2_cam1_Y)
pin <- pin[complete.cases(pin), ]
# get the number of pixels in the calibration
PixInPin <- (sqrt((pin$dp.pt1_cam1_X - pin$dp.pt2_cam1_X)^2 +
(pin$dp.pt1_cam1_Y-pin$dp.pt2_cam1_Y)^2)) /
metDat$CalSizeMM[ii] # to get to mm
# get anther and pollen locations
antherPoll <- data.frame(anthx = dp$pt3_cam1_X, anthy= dp$pt3_cam1_Y,
polx = dp$pt4_cam1_X, poly= dp$pt4_cam1_Y)
# get frame where pollen starts and leaves
antherPoll$polStart = 1:nrow(antherPoll) == metDat$framePollenStartedLeaving[ii]
antherPoll$polEnd = 1:nrow(antherPoll) == metDat$framePollenReleaseComplete[ii]
# gives only rows where either anth and pollen are complete
antherPoll = antherPoll[ complete.cases(antherPoll[c('anthx')]) |
complete.cases(antherPoll[c('polx')]), ]
# if x value starts to right of screen, reverse points,
# so all x values start on the left part of the screen at 0
if(lm(antherPoll[,1] ~ I(1:length( antherPoll[,1])))$coefficients[2] < 0 ){
antherPoll$anthx <- metDat[ii,'vidWidth'] - antherPoll$anthx
antherPoll$polx <- metDat[ii,'vidWidth'] - antherPoll$polx
}
# cbind data frame, to add smoothed columns
antherPoll <- data.frame(cbind(antherPoll, antherPoll))
# smooth with a smoothing spline, spar = 0.29 -- this smoothing
# spline was chosen through cross validation
foo = sapply(X = c("anthx.1", "anthy.1", "polx.1", "poly.1"), FUN = function(y){
yy = na.omit(antherPoll[, y])
xx = 1:length(yy)
sm1 <- smooth.spline(x = xx, y = yy, spar = 0.29)
antherPoll[, y][complete.cases(antherPoll[, y])] <<- sm1$y
})
# add time to data frame
antherPoll$tme = 0: (nrow(antherPoll) - 1) / fps # time
# add columns with absolute position into dataframe (calculated from smoothed data)
# calculate position from starting point, not from minimum point
bar = sapply(X = c("anthx.1", "anthy.1", "polx.1", "poly.1"), FUN = function(x){
newName = paste0(x, ".abs")
tmp <- antherPoll[,x] / PixInPin / 1000
antherPoll[,newName] <<- tmp - na.omit(tmp)[1]
#antherPoll[,newName] <<- tmp - min(na.omit(tmp))
})
# add columns to show velocity, based on smoothed, absolute position
# velocity is in m/s
bat = sapply(X = c("anthx.1.abs", "anthy.1.abs", "polx.1.abs", "poly.1.abs"),
FUN = function(x){
newName = paste0(x, ".vel")
tmp <- c(NaN, diff(antherPoll[,x])) * fps # add a NaN to beginning of data
antherPoll[,newName] <<- tmp
})
# calculate speed
antherPoll$anthspeed = sqrt(antherPoll$anthx.1.abs.vel^2 + antherPoll$anthy.1.abs.vel^2)
antherPoll$polspeed = sqrt(antherPoll$polx.1.abs.vel^2 + antherPoll$poly.1.abs.vel^2)
###########################################
# pollen acceleration
polVelocity = cbind(antherPoll$polx.1.abs.vel, antherPoll$poly.1.abs.vel)
polSpeed = antherPoll$polspeed
# plot(polSpeed)
tme = antherPoll$tme
polAccel = data.frame(rbind(c(NA, NA), apply(polVelocity, MARGIN = 2, FUN = diff))) * fps
# calculate unit tangent vector
T_t = polVelocity / polSpeed
DT = data.frame(rbind(c(NA, NA), apply(T_t, MARGIN = 2, FUN = diff))) * fps
NormDT = sqrt(DT[,1]^2 + DT[,2]^2)
Curvature = NormDT / polSpeed
# compute a_N (normal acceleration) and a_T (tangential acceleration)
# a_T = ds/dt
a_T = c(NA, diff(polSpeed) * fps)
N_t = data.frame(t(sapply(1:nrow(DT), FUN = function(x) unlist(DT[x, ] / NormDT[x]))))
# plot(a_T, type = "l", ylim = c(-3000, 3000))
# a_N = speed^2 * curvature
a_N = polSpeed^2 * Curvature
# check total accel by adding normal and tangential accelerations
# a_total = a_T * T_t + a_N * N_t
a_total = as.data.frame(t(sapply(X = 1:nrow(polAccel), FUN = function(x) a_T[x] * T_t[x, ] + a_N[x] * N_t[x,] )))
a_T_Pol = a_T
# calculate magnitude of acceleration, using two methods (should be the same)
# 1. Normal and tangential acceleration
a_mag1 = sqrt(a_T^2 + a_N^2)
amag2 = sqrt(polAccel[,1]^2 + polAccel[,2]^2)
###########################################
# anther acceleration
anthVelocity = cbind(antherPoll$anthx.1.abs.vel, antherPoll$anthy.1.abs.vel)
anthSpeed = antherPoll$anthspeed
tme = antherPoll$tme
anthAccel = data.frame(rbind(c(NA, NA), apply(anthVelocity, MARGIN = 2, FUN = diff))) * fps
# unit tangent vector
T_t = anthVelocity / anthSpeed
DT = data.frame(rbind(c(NA, NA), apply(T_t, MARGIN = 2, FUN = diff))) * fps
NormDT = sqrt(DT[,1]^2 + DT[,2]^2)
Curvature = NormDT / anthSpeed
# compute a_N (normal acceleration) and a_T (tangential acceleration)
# a_T = ds/dt
a_T = c(NA, diff(anthSpeed) * fps)
a_T_anth = a_T
N_t = data.frame(t(sapply(1:nrow(DT), FUN = function(x) unlist(DT[x, ] / NormDT[x]))))
# plot(a_T, type = "l", ylim = c(-3000, 3000))
# a_N = speed^2 * curvature
a_N = anthSpeed^2 * Curvature
# check total accel by adding normal and tangential accelerations
# a_total = a_T * T_t + a_N * N_t
a_total = as.data.frame(t(sapply(X = 1:nrow(anthAccel), FUN = function(x) a_T[x] * T_t[x, ] + a_N[x] * N_t[x,] )))
# align frames
tmeRoll <- seq(from = -which(antherPoll$polStart) + 1, length.out = length(tme)) / fps
# update for different visualization
# align frames -- center, based on max speed
# tmeRoll <- seq(from = -which.max(antherPoll$polspeed) + 1, length.out = length(tme)) / fps
dfi <- data.frame(anthSpeed, polSpeed, a_T_anth, a_T_Pol, tme,
trial = metDat$VideoName[ii],
tmeStart = antherPoll$polStart,
tmeEnd = antherPoll$polEnd,
centeredTime = tmeRoll)
newDF <- rbind(newDF,dfi)
print(ii)
}
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12
## [1] 13
## [1] 14
## [1] 15
## [1] 16
## [1] 17
## [1] 18
## [1] 19
## [1] 20
## [1] 21
## [1] 22
## [1] 23
## [1] 24
## [1] 25
## [1] 26
## [1] 27
## [1] 28
## [1] 29
## [1] 30
## [1] 31
## [1] 32
theme_set(theme_classic())
savePath = "/Users/callinswitzer/Dropbox/ExperSummer2016/Kalmia/Manuscript/FigsNotForSubmission/"
# anther speed
ggplot(newDF, aes(x = centeredTime, y = anthSpeed, group = trial)) +
geom_line(alpha = 0.5) +
xlim(c(-0.01, 0.02)) +
ylim(c(0,6)) +
labs(x = "Time (s)", y = "Anther speed (m/s)") +
ggsave(paste0(savePath, "antherSpeed01_CVSmoothSpline.pdf"), width = 5, height = 4)
## Warning: Removed 113 rows containing missing values (geom_path).
## Warning: Removed 113 rows containing missing values (geom_path).
# pollen speed
ggplot(newDF, aes(x = centeredTime, y = polSpeed, group = trial)) +
geom_line(alpha = 0.5) +
xlim(c(-0.01, 0.02)) +
ylim(c(0,6)) +
labs(x = "Time (s)", y = "Pollen speed (m/s)")
## Warning: Removed 53 rows containing missing values (geom_path).
ggsave(paste0(savePath, "pollenSpeed01_CVSmoothSpline.pdf"), width = 5, height = 4)
## Warning: Removed 53 rows containing missing values (geom_path).
# anther tangential acceleration
ggplot(newDF, aes(x = centeredTime, y = a_T_anth, group = trial)) +
geom_line(alpha = 0.5) +
#ylim(c(-2500, 4000)) +
xlim(c(-0.01, 0.02)) +
labs(x = "Time (s)", y = "Anther tangential acceleration (m/s/s)")
## Warning: Removed 145 rows containing missing values (geom_path).
ggsave(paste0(savePath, "antherTangAccel01_CVSmoothSpline.pdf"), width = 5, height = 4)
## Warning: Removed 145 rows containing missing values (geom_path).
# pollen tangential acceleration
# anther tangential acceleration
ggplot(newDF, aes(x = centeredTime, y = a_T_Pol, group = trial)) +
geom_line(alpha = 0.5) +
#ylim(c(-2500, 4000)) +
xlim(c(-0.01, 0.02)) +
labs(x = "Time (s)", y = "Pollen tangential acceleration (m/s/s)")
## Warning: Removed 85 rows containing missing values (geom_path).
ggsave(paste0(savePath, "PollenTangAccel01_CVSmoothSpline.pdf"), width = 5, height = 4)
## Warning: Removed 85 rows containing missing values (geom_path).
# anther speed
mmx <- as.data.frame(t(sapply(unique(as.character(newDF$trial)), FUN = function(x){
tmp <- newDF[newDF$trial == x, ]
return (unlist(tmp[which.max(tmp$anthSpeed),]))
})))
mmx$trial <- row.names(mmx)
ggplot() +
geom_line(data = newDF, aes(x = centeredTime, y = anthSpeed, group = as.factor(trial)), alpha = 0.5) +
xlim(c(-0.01, 0.02)) +
ylim(c(0,6)) +
labs(x = "Time (s)", y = "Anther speed (m/s)") +
geom_point(data = mmx, aes(x = centeredTime, y = anthSpeed), color = 'red', alpha = 0.5) +
theme(legend.position = "none")
## Warning: Removed 113 rows containing missing values (geom_path).
#+ facet_wrap(~ trial)
ggsave(paste0(savePath, "antherSpeedMax01_CVSmoothSpline.pdf"), width = 5, height = 4)
## Warning: Removed 113 rows containing missing values (geom_path).
# pollen speed
mmp <- as.data.frame(t(sapply(unique(as.character(newDF$trial)), FUN = function(x){
tmp <- newDF[newDF$trial == x, ]
tmp <- tmp[abs(tmp$centeredTime) < 0.01, ]
return (unlist(tmp[which.max(tmp$polSpeed),]))
})))
mmp$trial <- row.names(mmp)
# pollen speed
ggplot() +
geom_line(data = newDF, aes(x = centeredTime, y = polSpeed, group = trial), alpha = 0.5) +
xlim(c(-0.01, 0.02)) +
ylim(c(0,6)) +
labs(x = "Time (s)", y = "Pollen speed (m/s)") +
geom_point(data = mmp, aes(x = centeredTime, y = polSpeed), color = 'red', alpha = 0.5) +
theme(legend.position = "none")
## Warning: Removed 53 rows containing missing values (geom_path).
ggsave(paste0(savePath, "pollenSpeedMax01_CVSmoothSpline.pdf"), width = 5, height = 4)
## Warning: Removed 53 rows containing missing values (geom_path).
# anther acceleration
mma <- as.data.frame(t(sapply(unique(as.character(newDF$trial)), FUN = function(x){
tmp <- newDF[newDF$trial == x, ]
# get only points that are within 0.05 seconds of the centered time
# to ignore the anthers hitting the other side of the flower
#tmp <- tmp[abs(tmp$centeredTime) < 0.002, ]
return (unlist(tmp[which.max(tmp$a_T_anth),]))
})))
mma$trial <- row.names(mma)
ggplot() +
geom_line(data = newDF, aes(x = centeredTime, y = a_T_anth, group = trial), alpha = 0.5) +
#ylim(c(-2500, 4000)) +
xlim(c(-0.01, 0.02)) +
labs(x = "Time (s)", y = "Anther tangential acceleration (m/s/s)") +
geom_point(data = mma, aes(x = centeredTime, y = a_T_anth), color = 'red', alpha = 0.5)
## Warning: Removed 145 rows containing missing values (geom_path).
ggsave(paste0(savePath, "antherTangAccelMax01_CVSmoothSpline.pdf"), width = 5, height = 4)
## Warning: Removed 145 rows containing missing values (geom_path).
# pollen acceleration
mmpp <- as.data.frame(t(sapply(unique(as.character(newDF$trial)), FUN = function(x){
tmp <- newDF[newDF$trial == x, ]
#tmp <- tmp[abs(tmp$centeredTime) < 0.003, ]
return (unlist(tmp[which.max(tmp$a_T_Pol),]))
})))
mmpp$trial <- row.names(mmpp)
ggplot() +
geom_line(data = newDF, aes(x = centeredTime, y = a_T_Pol, group = trial), alpha = 0.5) +
#ylim(c(-2500, 4000)) +
xlim(c(-0.01, 0.02)) +
labs(x = "Time (s)", y = "Pollen tangential acceleration (m/s/s)") +
geom_point(data = mmpp, aes(x = centeredTime, y = a_T_Pol), color = 'red', alpha = 0.5)
## Warning: Removed 85 rows containing missing values (geom_path).
ggsave(paste0(savePath, "pollenTangAccelMax01_CVSmoothSpline.pdf"), width = 5, height = 4)
## Warning: Removed 85 rows containing missing values (geom_path).
ii = 28
print(ii)
## [1] 28
ddfile <- paste0(digdirect, metDat$redigitizedFile[ii])
ddfile <- "/Users/callinswitzer/Dropbox/ExperSummer2016/Kalmia/Manuscript/DatasetsSupplemental/CleanKalmiaDigitized/20150608_135726xypts.csv"
dp <- read.csv(ddfile)
# calibrate locations, based on digitized pin or other object
# calibration points
pin <- data.frame(dp$pt1_cam1_X, dp$pt1_cam1_Y, dp$pt2_cam1_X, dp$pt2_cam1_Y)
pin <- pin[complete.cases(pin), ]
# get the number of pixels in the calibration
PixInPin <- (sqrt((pin$dp.pt1_cam1_X - pin$dp.pt2_cam1_X)^2 +
(pin$dp.pt1_cam1_Y-pin$dp.pt2_cam1_Y)^2)) /
metDat$CalSizeMM[ii] # to get to mm
# get anther and pollen locations
antherPoll <- data.frame(anthx = dp$pt3_cam1_X, anthy= dp$pt3_cam1_Y,
polx = dp$pt4_cam1_X, poly= dp$pt4_cam1_Y)
# get frame where pollen starts and leaves
antherPoll$polStart = 1:nrow(antherPoll) == metDat$framePollenStartedLeaving[ii]
antherPoll$polEnd = 1:nrow(antherPoll) == metDat$framePollenReleaseComplete[ii]
# gives only rows where either anth and pollen are complete
antherPoll = antherPoll[ complete.cases(antherPoll[c('anthx')]) |
complete.cases(antherPoll[c('polx')]), ]
# if x value starts to right of screen, reverse points,
# so all x values start on the left part of the screen at 0
if(lm(antherPoll[,1] ~ I(1:length( antherPoll[,1])))$coefficients[2] < 0 ){
antherPoll$anthx <- metDat[ii,'vidWidth'] - antherPoll$anthx
antherPoll$polx <- metDat[ii,'vidWidth'] - antherPoll$polx
}
plot(antherPoll$polx, antherPoll$poly, asp = 1)
# cbind data frame, to add smoothed columns
antherPoll <- data.frame(cbind(antherPoll, antherPoll))
# smooth with a smoothing spline, spar = 0.29 -- this smoothing
# spline was chosen through cross validation
foo = sapply(X = c("anthx.1", "anthy.1", "polx.1", "poly.1"), FUN = function(y){
yy = na.omit(antherPoll[, y])
xx = 1:length(yy)
sm1 <- smooth.spline(x = xx, y = yy, spar = 0.29)
antherPoll[, y][complete.cases(antherPoll[, y])] <<- sm1$y
})
plot(antherPoll$anthx.1, antherPoll$anthy.1)
antherPoll <- antherPoll[10:nrow(antherPoll), ]
# add time to data frame
antherPoll$tme = 0: (nrow(antherPoll) - 1) / fps # time
# add columns with absolute position into dataframe (calculated from smoothed data)
# calculate position from starting point, not from minimum point
bar = sapply(X = c("anthx.1", "anthy.1", "polx.1", "poly.1"), FUN = function(x){
newName = paste0(x, ".abs")
tmp <- antherPoll[,x] / PixInPin / 1000
antherPoll[,newName] <<- tmp - na.omit(tmp)[1]
#antherPoll[,newName] <<- tmp - min(na.omit(tmp))
})
plot(antherPoll$anthx.1.abs, antherPoll$anthy.1.abs)
# add columns to show velocity, based on smoothed, absolute position
# velocity is in m/s
bat = sapply(X = c("anthx.1.abs", "anthy.1.abs", "polx.1.abs", "poly.1.abs"),
FUN = function(x){
newName = paste0(x, ".vel")
tmp <- c(NaN, diff(antherPoll[,x])) * fps # add a NaN to beginning of data
antherPoll[,newName] <<- tmp
})
# calculate speed
antherPoll$anthspeed = sqrt(antherPoll$anthx.1.abs.vel^2 + antherPoll$anthy.1.abs.vel^2)
antherPoll$polspeed = sqrt(antherPoll$polx.1.abs.vel^2 + antherPoll$poly.1.abs.vel^2)
###########################################
# pollen acceleration
polVelocity = cbind(antherPoll$polx.1.abs.vel, antherPoll$poly.1.abs.vel)
polSpeed = antherPoll$polspeed
# plot(polSpeed)
tme = antherPoll$tme
polAccel = data.frame(rbind(c(NA, NA), apply(polVelocity, MARGIN = 2, FUN = diff))) * fps
# calculate unit tangent vector
T_t = polVelocity / polSpeed
DT = data.frame(rbind(c(NA, NA), apply(T_t, MARGIN = 2, FUN = diff))) * fps
NormDT = sqrt(DT[,1]^2 + DT[,2]^2)
Curvature = NormDT / polSpeed
# compute a_N (normal acceleration) and a_T (tangential acceleration)
# a_T = ds/dt
a_T = c(NA, diff(polSpeed) * fps)
N_t = data.frame(t(sapply(1:nrow(DT), FUN = function(x) unlist(DT[x, ] / NormDT[x]))))
# plot(a_T, type = "l", ylim = c(-3000, 3000))
# a_N = speed^2 * curvature
a_N = polSpeed^2 * Curvature
# check total accel by adding normal and tangential accelerations
# a_total = a_T * T_t + a_N * N_t
a_total = as.data.frame(t(sapply(X = 1:nrow(polAccel), FUN = function(x) a_T[x] * T_t[x, ] + a_N[x] * N_t[x,] )))
a_T_Pol = a_T
ggplot(antherPoll, aes(y = poly.1.abs, x = polx.1.abs)) +
geom_path(alpha = 0.5) +
coord_fixed() +
ylim(0, 0.11)+
labs(y = "Vertical Position (m)", x = "Horizontal position (m)")
x1 <- ggplot(antherPoll, aes(x = tme, y = polx.1.abs)) +
geom_line() +
labs(x = "", y = "Pollen\nhorizontal position (m)") +
ylim(-0.01, 0.11)+
theme(axis.title.y=element_text(margin=margin(0,20,0,0)))+
theme(plot.margin = unit(c(1,1,-1,0), "cm")) +
ggsave(filename = paste0(savePath,"/horizontalpos.png"), width = 5*1.5, height = 1.5*1.5, dpi =200)
x2 <- ggplot(antherPoll, aes(x = tme, y = poly.1.abs)) +
geom_line() +
ylim(-0.01, 0.11)+
labs(x = "", y = "Pollen\nvertical Position (m)")+
theme(axis.title.y=element_text(margin=margin(0,20,0,0)))+
theme(plot.margin = unit(c(1,1,-1,0), "cm")) +
ggsave(filename = paste0(savePath,"/verticalpos.png"), width = 5*1.5, height = 1.5*1.5, dpi =200)
x3 <- ggplot(antherPoll, aes(x = tme, y = polspeed)) +
geom_line() +
labs(x = "", y = "Pollen\nspeed (m/s)")+
theme(plot.margin = unit(c(1,1,-5,0), "cm")) +
theme(axis.title.y=element_text(margin=margin(0,20,0,0)))
x3
## Warning: Removed 1 rows containing missing values (geom_path).
ggsave(filename = paste0(savePath,"/speed.png"), width = 5*1.5, height = 1.5*1.5, dpi =200)
## Warning: Removed 1 rows containing missing values (geom_path).
x4 <- ggplot(antherPoll, aes(x = tme, y = a_T_Pol)) +
geom_line() +
labs(x = "Time (s)", y = "Pollen\ntangential accel. (m/s/s)") +
theme(axis.title.y=element_text(margin=margin(0,20,0,0)))+
theme(plot.margin = unit(c(0,0,0,0), "cm")) +
ggsave(filename = paste0(savePath,"/pollenAccel.png"), width = 5*1.5, height = 1.5*1.5, dpi =200)
## Warning: Removed 2 rows containing missing values (geom_path).
#install.packages("cowplot")
library(cowplot)
##
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggplot2':
##
## ggsave
plot_grid(x1, x2, x3, x4, align='vh', ncol = 1)
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 2 rows containing missing values (geom_path).
ggsave(paste0(savePath,"fourplots.png"), height = 10, width = 10, dpi = 120)
theme_set(theme_classic())
savePath = "/Users/callinswitzer/Dropbox/ExperSummer2016/Kalmia/Manuscript/FigsNotForSubmission/"
# anther speed
ggplot(newDF, aes(x = centeredTime, y = anthSpeed, group = trial)) +
geom_line(alpha = 0.5) +
xlim(c(-0.01, 0.02)) +
ylim(c(0,6)) +
labs(x = "Time (s)", y = "Anther speed (m/s)") +
ggsave(paste0(savePath, "antherSpeed01_CVSmoothSpline.pdf"), width = 5, height = 4)
## Warning: Removed 113 rows containing missing values (geom_path).
## Warning: Removed 113 rows containing missing values (geom_path).
# pollen speed
ggplot(newDF, aes(x = centeredTime, y = polSpeed, group = trial)) +
geom_line(alpha = 0.5) +
xlim(c(-0.01, 0.02)) +
ylim(c(0,6)) +
labs(x = "Time (s)", y = "Pollen speed (m/s)")
## Warning: Removed 53 rows containing missing values (geom_path).
ggsave(paste0(savePath, "pollenSpeed01_CVSmoothSpline.pdf"), width = 5, height = 4)
## Warning: Removed 53 rows containing missing values (geom_path).
# anther tangential acceleration
ggplot(newDF, aes(x = centeredTime, y = a_T_anth, group = trial)) +
geom_line(alpha = 0.5) +
#ylim(c(-2500, 4000)) +
xlim(c(-0.01, 0.02)) +
labs(x = "Time (s)", y = "Anther tangential acceleration (m/s/s)")
## Warning: Removed 145 rows containing missing values (geom_path).
ggsave(paste0(savePath, "antherTangAccel01_CVSmoothSpline.pdf"), width = 5, height = 4)
## Warning: Removed 145 rows containing missing values (geom_path).
# pollen tangential acceleration
# anther tangential acceleration
ggplot(newDF, aes(x = centeredTime, y = a_T_Pol, group = trial)) +
geom_line(alpha = 0.5) +
#ylim(c(-2500, 4000)) +
xlim(c(-0.01, 0.02)) +
labs(x = "Time (s)", y = "Pollen tangential acceleration (m/s/s)")
## Warning: Removed 85 rows containing missing values (geom_path).
ggsave(paste0(savePath, "PollenTangAccel01_CVSmoothSpline.pdf"), width = 5, height = 4)
## Warning: Removed 85 rows containing missing values (geom_path).
# anther speed
mmx <- as.data.frame(t(sapply(unique(as.character(newDF$trial)), FUN = function(x){
tmp <- newDF[newDF$trial == x, ]
return (unlist(tmp[which.max(tmp$anthSpeed),]))
})))
mmx$trial <- row.names(mmx)
ggplot() +
geom_line(data = newDF, aes(x = centeredTime, y = anthSpeed, group = as.factor(trial)), alpha = 0.5) +
xlim(c(-0.01, 0.02)) +
ylim(c(0,6)) +
labs(x = "Time (s)", y = "Anther speed (m/s)") +
geom_point(data = mmx, aes(x = centeredTime, y = anthSpeed), color = 'red', alpha = 0.5) +
theme(legend.position = "none")
## Warning: Removed 113 rows containing missing values (geom_path).
#+ facet_wrap(~ trial)
ggsave(paste0(savePath, "antherSpeedMax01_CVSmoothSpline.pdf"), width = 5, height = 4)
## Warning: Removed 113 rows containing missing values (geom_path).
# pollen speed
mmp <- as.data.frame(t(sapply(unique(as.character(newDF$trial)), FUN = function(x){
tmp <- newDF[newDF$trial == x, ]
tmp <- tmp[abs(tmp$centeredTime) < 0.01, ]
return (unlist(tmp[which.max(tmp$polSpeed),]))
})))
mmp$trial <- row.names(mmp)
md = merge(x = mmx[, c('trial', 'anthSpeed')], metDat, by.x = "trial", by.y =
"VideoName")
md = merge(x = mmp[, c('trial', 'polSpeed')], md, by = "trial")
md = merge(x = mmpp[, c('trial', 'a_T_Pol')], md, by = "trial")
md = merge(x = mma[, c('trial', 'a_T_anth')], md, by = "trial")
#LMER
hist(md$anthSpeed) #dist is fine
modVelMaxAnth <- lmer(formula = anthSpeed ~ (1|plant/FlowerNumber), data = md)
summary(modVelMaxAnth)
## Linear mixed model fit by REML ['lmerMod']
## Formula: anthSpeed ~ (1 | plant/FlowerNumber)
## Data: md
##
## REML criterion at convergence: 63.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.52444 -0.65641 0.00419 0.52344 1.89172
##
## Random effects:
## Groups Name Variance Std.Dev.
## FlowerNumber:plant (Intercept) 0.2844 0.5333
## plant (Intercept) 0.4103 0.6406
## Residual 0.1608 0.4010
## Number of obs: 32, groups: FlowerNumber:plant, 15; plant, 7
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.772 0.296 12.74
confint(modVelMaxAnth)
## Computing profile confidence intervals ...
## 2.5 % 97.5 %
## .sig01 0.2798029 0.9608723
## .sig02 0.0000000 1.2573804
## .sigma 0.2973171 0.5798682
## (Intercept) 3.1531345 4.3868581
# diagnostics
plot(modVelMaxAnth)
# QQPlot for group-level effects
qqnorm(ranef(modVelMaxAnth)$plant[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(modVelMaxAnth)$plant[[1]])
# QQPlot for group-level effects
qqnorm(ranef(modVelMaxAnth)$FlowerNumber[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(modVelMaxAnth)$FlowerNumber[[1]]) #
modVelMaxPol <- lmer(formula = polSpeed ~ (1|plant/FlowerNumber), data = md)
summary(modVelMaxPol) # final model for paper
## Linear mixed model fit by REML ['lmerMod']
## Formula: polSpeed ~ (1 | plant/FlowerNumber)
## Data: md
##
## REML criterion at convergence: 71.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.03862 -0.49111 -0.07546 0.27208 1.95680
##
## Random effects:
## Groups Name Variance Std.Dev.
## FlowerNumber:plant (Intercept) 0.2414 0.4914
## plant (Intercept) 0.1210 0.3479
## Residual 0.3089 0.5558
## Number of obs: 32, groups: FlowerNumber:plant, 15; plant, 7
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.5019 0.2141 16.36
confint(modVelMaxPol)
## Computing profile confidence intervals ...
## 2.5 % 97.5 %
## .sig01 0.06228879 0.9128825
## .sig02 0.00000000 0.8681086
## .sigma 0.41375779 0.7940168
## (Intercept) 3.05487009 3.9520846
plot(modVelMaxPol)
# log transformed, b/c distribution is skewed.
hist(md$a_T_Pol)
hist(log(md$a_T_Pol)) #better
modAccMaxPol <- lmer(formula = log(a_T_Pol) ~ (1|plant/FlowerNumber), data = md)
summary(modAccMaxPol) #final model for paper
## Linear mixed model fit by REML ['lmerMod']
## Formula: log(a_T_Pol) ~ (1 | plant/FlowerNumber)
## Data: md
##
## REML criterion at convergence: 13.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.62436 -0.66206 0.03069 0.61934 2.08762
##
## Random effects:
## Groups Name Variance Std.Dev.
## FlowerNumber:plant (Intercept) 0.001824 0.0427
## plant (Intercept) 0.066778 0.2584
## Residual 0.056679 0.2381
## Number of obs: 32, groups: FlowerNumber:plant, 15; plant, 7
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 8.3302 0.1095 76.1
exp(confint(modAccMaxPol)) # CI for paper
## Computing profile confidence intervals ...
## 2.5 % 97.5 %
## .sig01 1.000000 1.310279
## .sig02 1.000000 1.646929
## .sigma 1.194421 1.392101
## (Intercept) 3308.135397 5250.567110
plot(modAccMaxPol)
# QQPlot for group-level effects
qqnorm(ranef(modAccMaxPol)$plant[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(modAccMaxPol)$plant[[1]])
# QQPlot for group-level effects
qqnorm(ranef(modAccMaxPol)$FlowerNumber[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(modAccMaxPol)$FlowerNumber[[1]])
# not reporting anthers, because they're basically the same as the pollen
modAccMaxAnth <- lmer(formula = log(a_T_anth) ~ (1|plant/FlowerNumber), data = md)
summary(modAccMaxAnth)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log(a_T_anth) ~ (1 | plant/FlowerNumber)
## Data: md
##
## REML criterion at convergence: 19
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.57600 -0.38059 0.05996 0.30044 2.00083
##
## Random effects:
## Groups Name Variance Std.Dev.
## FlowerNumber:plant (Intercept) 0.01712 0.1308
## plant (Intercept) 0.10281 0.3206
## Residual 0.05554 0.2357
## Number of obs: 32, groups: FlowerNumber:plant, 15; plant, 7
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 8.2824 0.1355 61.12
exp(confint(modAccMaxAnth))
## Computing profile confidence intervals ...
## 2.5 % 97.5 %
## .sig01 1.000000 1.408217
## .sig02 1.024878 1.843027
## .sigma 1.190678 1.405926
## (Intercept) 2988.212500 5284.547669
plot(modAccMaxAnth)
# diagnostics
plot(modAccMaxAnth)
# QQPlot for group-level effects
qqnorm(ranef(modAccMaxAnth)$plant[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(modAccMaxAnth)$plant[[1]])
# QQPlot for group-level effects
qqnorm(ranef(modAccMaxAnth)$FlowerNumber[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(modAccMaxAnth)$FlowerNumber[[1]]) #